sampleQC <- function(bamFile,bedFile=NULL,blklist=NULL,ChrOfInterest=NULL,GeneAnnotation=NULL,Window=400,FragmentLength=50,shiftWindowStart=1,shiftWindowEnd=2,mapQCutoff=15){
#    require(Rsamtools)
#    require(GenomicRanges)
#    require(GenomicAlignments)
   
   ChrLengths <- scanBamHeader(bamFile)[[1]]$targets
   message("Bam file has ",length(names(ChrLengths))," contigs")
   if(!all(ChrOfInterest %in% names(ChrLengths))){
      stop("Contigs of interest are not all in Bam file!!")
   }
   if(!is.null(ChrOfInterest)){
      ChrLengths <- ChrLengths[names(ChrLengths) %in% ChrOfInterest]
   }
   
   ShiftMat <- NULL
   FlagTagCounts <- NULL
   CovHist <- NULL
   Cov <- NULL
   SSD <- NULL
   readlength <- as.numeric(NA)
   CoverageMatrix <- NULL
   CoverageMatrixInput <- NULL
   bedRangesTemp <- GetGRanges(GRanges(),AllChr=names(ChrLengths))
   bedRangesSummitsTemp <- vector("numeric")
   Counts <- NULL
   MultiRangesCountsTemp <- NULL
   txdb <- NULL
   GeneAnnotationTotalSizes <- NULL
   if(class(GeneAnnotation)!="list" & !is.null(GeneAnnotation)) {
      GeneAnnotation = getAnnotation(GeneAnnotation,AllChr=names(ChrLengths))
 
   }
   if(class(GeneAnnotation)=="list"){
      GSizes <- unlist(lapply(GeneAnnotation,function(x)sum(width(x))))
      GeneAnnotationTotalSizes <- GSizes[!names(GSizes) %in% "version"] 
   }
   
   for(k in 1:length(ChrLengths)){
      
      Param <- ScanBamParam(which=GRanges(seqnames=names(ChrLengths)[k],IRanges(start=1,end=unname(ChrLengths[names(ChrLengths) == names(ChrLengths)[k]])-shiftWindowEnd)),
                            what=c("flag","mapq"))
      temp <- readGAlignmentsFromBam(bamFile,param=Param)
      if(length(temp) > 1){
      if(k == 1){
        tocheckforreads <- 1000
        readlength=round(mean(width(temp[1:tocheckforreads])))
      }
      
      Sample_GIT <- GIntervalTree(GRanges(seqnames=seqnames(temp),ranges=ranges(temp),strand=strand(temp),elementMetadata(temp)))
      
      flagMapQ <- cbind(bamFlagAsBitMatrix(Sample_GIT$flag)[,c("isUnmappedQuery","isDuplicate")],Sample_GIT$mapq)
      colnames(flagMapQ) <- c("A","B","C")
      temp <- as.data.frame(xtabs(~A+B+C, data=flagMapQ))
      UnMapped <- sum(temp[temp[,"A"] == 1,"Freq"])
      Mapped <- sum(temp[temp[,"A"] != 1,"Freq"])
      Duplicates <- sum(temp[temp[,"A"] != 1 & temp[,"B"] == 1,"Freq"])
      
      MapQPass <- sum(temp[temp[,"A"] != 1 & as.numeric(temp[,"C"]) >= 15,"Freq"])
      MapQPassAndDup <- sum(temp[temp[,"A"] != 1 & temp[,"B"] == 1 & as.numeric(temp[,"C"]) >= mapQCutoff,"Freq"])
      FlagTagCounts <- rbind(cbind(UnMapped,Mapped,Duplicates,MapQPass,MapQPassAndDup),FlagTagCounts)
      
      
      Cov <- coverage(Sample_GIT)
      message("Calculating coverage histogram for ",names(ChrLengths)[k],"\n")
      CovHist <- c(CovHist,list(colSums(table(Cov))))
      message("Calculating SSD for ",names(ChrLengths)[k],"\n")
      SSD <- c(SSD,sd(Cov))
      
      PosCoverage <- coverage(IRanges(start(Sample_GIT[strand(Sample_GIT)=="+"]),start(Sample_GIT[strand(Sample_GIT)=="+"])),width=ChrLengths[k])
      NegCoverage <- coverage(IRanges(end(Sample_GIT[strand(Sample_GIT)=="-"]),end(Sample_GIT[strand(Sample_GIT)=="-"])),width=ChrLengths[k])
      message("Calculating shift for ",names(ChrLengths)[k],"\n")
      #ShiftsTemp <- shiftApply(seq(shiftWindowStart,shiftWindowEnd),PosCoverage,NegCoverage,cor)
      ShiftsTemp <- shiftApply(seq(shiftWindowStart,shiftWindowEnd),PosCoverage,NegCoverage,RleSumAny, verbose = T)
      ShiftMat <- cbind(ShiftMat,ShiftsTemp)
      
      if(!is.null(bedFile)){
        bedRanges <- GetGRanges(bedFile,as.vector(names(ChrLengths)),names(ChrLengths)[k])
        GRangesOfInterestList <- GRangesList(GRanges(seqnames(bedRanges),ranges(bedRanges)))
        names(GRangesOfInterestList) <- c(names(GRangesOfInterestList),"Peaks")
      }else{
        GRangesOfInterestList <- GRangesList()
      }
      
      if(!is.null(blklist)){
         blkRanges <- GetGRanges(blklist,names(ChrLengths),names(ChrLengths)[k])           
         GRangesOfInterestListBL <- GRangesList(GRanges(seqnames(blkRanges),ranges(blkRanges)))
         names(GRangesOfInterestListBL) <- "BlackList"
         GRangesOfInterestList <- c(GRangesOfInterestList,GRangesOfInterestListBL)
         #print(names(GRangesOfInterestList))
      }
      if(class(GeneAnnotation)=="list"){

         GRangesOfInterestListGA <- 
                                    GRangesList(GetGRanges(GeneAnnotation$All5utrs,names(ChrLengths),names(ChrLengths)[k],simplify=T),
                                                GetGRanges(GeneAnnotation$All3utrs,names(ChrLengths),names(ChrLengths)[k],simplify=T),
                                                GetGRanges(GeneAnnotation$Allcds,names(ChrLengths),names(ChrLengths)[k],simplify=T),
                                                GetGRanges(GeneAnnotation$Allintrons,names(ChrLengths),names(ChrLengths)[k],simplify=T),
                                                GetGRanges(GeneAnnotation$Alltranscripts,names(ChrLengths),names(ChrLengths)[k],simplify=T),
                                                GetGRanges(GeneAnnotation$Promoters500,names(ChrLengths),names(ChrLengths)[k],simplify=T),
                                                GetGRanges(GeneAnnotation$Promoters2000to500,names(ChrLengths),names(ChrLengths)[k],simplify=T),
                                                GetGRanges(GeneAnnotation$LongPromoter20000to2000,names(ChrLengths),names(ChrLengths)[k],simplify=T))
                                                names(GRangesOfInterestListGA) <-
                                                c("All5utrs","All3utrs","Allcds","Allintrons","Alltranscripts",
                                                "Promoters500","Promoters2000to500","LongPromoter20000to2000")
          
         GRangesOfInterestList <- c(GRangesOfInterestList,GRangesOfInterestListGA)
      }
      MultiRangesCountsTemp <- c(MultiRangesCountsTemp,countOverlaps(GRangesOfInterestList,Sample_GIT))
      message("Counting reads in features for ",names(ChrLengths)[k],"\n")        
      if(!is.null(bedFile)){
        CountsTemp <- countOverlaps(bedRanges,Sample_GIT)
        Counts  <- c(Counts,CountsTemp)
        
        bedRangesTemp <- c(bedRangesTemp,bedRanges)
        message("Signal over peaks for ",names(ChrLengths)[k],"\n")                   
        AllFragRanges <- resize(as(Sample_GIT[Sample_GIT %over% bedRanges],"GRanges"),FragmentLength,"start")
        bedRangesSummits <- findCovMaxPos(AllFragRanges,bedRanges,names(ChrLengths)[k],FragmentLength)
        bedRangesSummitsTemp <- c(bedRangesSummitsTemp,as.numeric(as.vector(start(bedRangesSummits))))
        updatedRanges  <- resize(bedRangesSummits,Window,"center")
        AllCov <- coverage(AllFragRanges,width=unname(ChrLengths[k])+Window)
        CoverageMatrix <- rbind(CoverageMatrix,matrix(as.vector(AllCov[[which(names(AllCov) == names(ChrLengths)[k])]][ranges(updatedRanges[seqnames(updatedRanges) == names(ChrLengths)[k]])]),ncol=Window,byrow=T))
      }
   }
   }
   if(!is.null(FlagTagCounts)){
    FlagTagCounts <- colSums(FlagTagCounts)
   }else{
     FlagTagCounts <- as.numeric(data.frame(UnMapped=0,Mapped=0,Duplicates=0,MapQPass=0,MapQPassAndDup=0))
   }
   Weights <- ChrLengths
   CovHistAll <- NULL
   if(!is.null(CovHist)){
   for(l in 1:length(CovHist)){
       print(length(l))
       tempCovHist <- cbind(as.numeric(names(CovHist[[l]])),as.numeric(CovHist[[l]]))
       tempCovHist <- merge(CovHistAll,tempCovHist,by.x=0,by.y=1,all=T,sort=F)
       CovSums <- data.frame(Depth=rowSums(tempCovHist[,-1,drop=F],na.rm=T),row.names=tempCovHist[,1])
       CovHistAll <- CovSums
     }
     tempCovHistAll <- as.numeric(CovHistAll[,1])
     names(tempCovHistAll) <- rownames(CovHistAll)
     CovHistAll <- tempCovHistAll
   }else{
     CovHistAll <- as.numeric(NA)     
   }
   if(!is.null(ShiftMat)){
   ShiftsAv <- apply(ShiftMat,1,function(x)weighted.mean(x,Weights,na.rm=T))
   }else{
     ShiftsAv <- as.numeric(NA)
   }
   if(!is.null(SSD)){   
   SSDAv <- (weighted.mean(sd(Cov)[names(ChrLengths)],ChrLengths)*1000)/sqrt(FlagTagCounts[4])
  }else{
    SSDAv <- as.numeric(NA)     
  }

   if(!is.null(MultiRangesCountsTemp)){
     GFCountsMatrix <- matrix(MultiRangesCountsTemp,ncol=length(GRangesOfInterestList),byrow=T)
     colnames(GFCountsMatrix) <- unique(names(MultiRangesCountsTemp))
   }else{
     GFCountsMatrix <- NULL
   }
   
   if(!is.null(bedFile)){
     AvProfile <- colMeans(CoverageMatrix)
     NormAvProfile <- (AvProfile/FlagTagCounts[4])*1e6
     elementMetadata(bedRangesTemp) <- data.frame(Counts,bedRangesSummitsTemp)
     #print(length(GRangesOfInterestList))
     ReadsInPeaks <- sum(GFCountsMatrix[,"Peaks"])
   }else{
     AvProfile <- NA
     NormAvProfile <- NA
     bedRangesTemp <- GRanges()
     #print(length(GRangesOfInterestList))
     ReadsInPeaks <- as.numeric(NA)
   }  
   if(!is.null(blklist)){   
    ReadsInBLs <- sum(GFCountsMatrix[,"BlackList"])  
   }else{
    ReadsInBLs <- as.numeric(NA)
   }

   if(class(GeneAnnotation)=="list" & !is.null(GFCountsMatrix)){
     ReadsIn5utrs<- sum(GFCountsMatrix[,"All5utrs"])
     ReadsIn3utrs <- sum(GFCountsMatrix[,"All3utrs"])
     ReadsInCDS <- sum(GFCountsMatrix[,"Allcds"])
     ReadsInIntrons <- sum(GFCountsMatrix[,"Allintrons"])
     ReadsInTranscripts <- sum(GFCountsMatrix[,"Alltranscripts"])              
     ReadsInPromoters500 <- sum(GFCountsMatrix[,"Promoters500"])              
     ReadsInPromoters2000to500 <- sum(GFCountsMatrix[,"Promoters2000to500"])              
     ReadsInLongPromoter20000to2000 <- sum(GFCountsMatrix[,"LongPromoter20000to2000"])              
     CountsInFeatures <- list(CountsIn5utrs=ReadsIn5utrs,CountsIn3UTRs=ReadsIn3utrs,
     CountsinReads=ReadsInCDS,CountsInIntrons=ReadsInIntrons,CountsInTranscripts=ReadsInTranscripts,
     CountsInPromoters500=ReadsInPromoters500,CountsInPromoters2000to500=ReadsInPromoters2000to500,
     CountsInLongPromoter20000to2000=ReadsInLongPromoter20000to2000)

     TotalLength5utrs<- GeneAnnotationTotalSizes["All5utrs"]
     TotalLength3utrs <- GeneAnnotationTotalSizes["All3utrs"]
     TotalLengthCDS <- GeneAnnotationTotalSizes["Allcds"]
     TotalLengthIntrons <- GeneAnnotationTotalSizes["Allintrons"]
     TotalLengthTranscripts <- GeneAnnotationTotalSizes["Alltranscripts"]              
     TotalLengthPromoters500 <- GeneAnnotationTotalSizes["Promoters500"]              
     TotalLengthPromoters2000to500 <- GeneAnnotationTotalSizes["Promoters2000to500"]              
     TotalLengthLongPromoter20000to2000 <- GeneAnnotationTotalSizes["LongPromoter20000to2000"]
     PropInFeatures <- list(PropIn5utrs=TotalLength5utrs/sum(as.numeric(ChrLengths)),PropIn3UTRs=TotalLength3utrs/sum(as.numeric(ChrLengths)),
                              PropInCDS=TotalLengthCDS/sum(as.numeric(ChrLengths)),PropInIntrons=TotalLengthIntrons/sum(as.numeric(ChrLengths)),PropInTranscripts=TotalLengthTranscripts/sum(as.numeric(ChrLengths)),
                              PropInPromoters500=TotalLengthPromoters500/sum(as.numeric(ChrLengths)),PropInPromoters2000to500=TotalLengthPromoters2000to500/sum(as.numeric(ChrLengths)),
                              PropInLongPromoter20000to2000=TotalLengthLongPromoter20000to2000/sum(as.numeric(ChrLengths))
    )
     
     
     
     
   }else{
     ReadsIn5utrs<- NA
     ReadsIn3utrs <- NA
     ReadsInCDS <- NA
     ReadsInIntrons <- NA
     ReadsInTranscripts <- NA
     ReadsInPromoters500 <- NA
     ReadsInPromoters2000to500 <- NA
     ReadsInLongPromoter20000to2000 <- NA

     CountsInFeatures <- list(CountsIn5utrs=ReadsIn5utrs,CountsIn3UTRs=ReadsIn3utrs,
     CountsinReads=ReadsInCDS,CountsInIntrons=ReadsInIntrons,CountsInTranscripts=ReadsInTranscripts,
     CountsInPromoters500=ReadsInPromoters500,CountsInPromoters2000to500=ReadsInPromoters2000to500,
     CountsInLongPromoter20000to2000=ReadsInLongPromoter20000to2000) 

     
     TotalLength5utrs<- NA
     TotalLength3utrs <- NA
     TotalLengthCDS <- NA
     TotalLengthIntrons <- NA
     TotalLengthTranscripts <- NA
     TotalLengthPromoters500 <- NA
     TotalLengthPromoters2000to500 <- NA
     TotalLengthLongPromoter20000to2000 <- GeneAnnotationTotalSizes["LongPromoter20000to2000"]
     PropInFeatures <- list(PropIn5utrs=TotalLength5utrs,PropIn3UTRs=TotalLength3utrs,
                            PropInReads=TotalLengthCDS,PropInIntrons=TotalLengthIntrons,PropInTranscripts=TotalLengthTranscripts,
                            PropInPromoters500=TotalLengthPromoters500,PropInPromoters2000to500=TotalLengthPromoters2000to500,
                            PropInLongPromoter20000to2000=TotalLengthLongPromoter20000to2000)
     
     
   }
   
   
   
   CHQC <- new("ChIPQCsample", bedRangesTemp,AveragePeakSignal=list(AvProfile,NormAvProfile),
               CrossCoverage=ShiftsAv,SSD=SSDAv,CountsInPeaks=ReadsInPeaks,
               CountsInBlackList=ReadsInBLs,
               CountsInFeatures=CountsInFeatures,PropInFeatures=PropInFeatures,
               CoverageHistogram=CovHistAll,FlagAndTagCounts=FlagTagCounts,readlength=readlength)

   return(CHQC)    
   
}

RleSumAny <- function (e1, e2)
{
   len <- length(e1)
   stopifnot(len == length(e2))
   x1 <- runValue(e1); s1 <- cumsum(runLength(e1))
   x2 <- runValue(e2); s2 <- cumsum(runLength(e2))
   .Call("rle_sum_any",
         as.integer(x1), as.integer(s1),
         as.integer(x2), as.integer(s2),
         as.integer(len),
         PACKAGE = "chipseq")
}

GetGRanges <- function(LoadFile,AllChr=NULL,ChrOfInterest=NULL,simple=F,sepr="\t",simplify=F){
#    require(Rsamtools)
#    require(GenomicRanges)
   
   if(class(LoadFile) == "GRanges"){
      RegionRanges <- LoadFile
      if(simplify){
         RegionRanges <- GRanges(seqnames(RegionRanges),ranges(RegionRanges))
      }
   }else{
      if(class(LoadFile) == "character"){
         RangesTable <- read.delim(LoadFile,sep=sepr,header=T,skip="#")
      }else if(class(LoadFile) == "matrix"){
         RangesTable <- as.data.frame(LoadFile)
      } else{
         RangesTable <- as.data.frame(LoadFile)
      }
      Chromosomes <- as.vector(RangesTable[,1])
      Start <- as.numeric(as.vector(RangesTable[,2]))
      End <- as.numeric(as.vector(RangesTable[,3]))
      RegionRanges <- GRanges(seqnames=Chromosomes,ranges=IRanges(start=Start,end=End))
      if(simple == F){
         if(ncol(RangesTable) > 4){
            ID <- as.vector(RangesTable[,4])
            Score <- as.vector(RangesTable[,5])
            if(ncol(RangesTable) > 6){
               Strand <- rep("*",nrow(RangesTable))
               RemainderColumn <- as.data.frame(RangesTable[,-c(1:6)])
               elementMetadata(RegionRanges) <- cbind(ID,Score,Strand,RemainderColumn)
            }else{
               elementMetadata(RegionRanges) <- cbind(ID,Score)
            }
         }
      }
   }
   if(!is.null(AllChr)){ 
      RegionRanges <- RegionRanges[seqnames(RegionRanges) %in% AllChr]    
      seqlevels(RegionRanges,force=T) <- AllChr
   }
   if(!is.null(ChrOfInterest)){      
      RegionRanges <- RegionRanges[seqnames(RegionRanges) == ChrOfInterest]      
   }
   
   return(RegionRanges)
}

findCovMaxPos <- function(reads,bedRanges,ChrOfInterest,FragmentLength){
#    require(GenomicRanges)
#    require(Rsamtools)
   
   cat("done\n")
   cat("Calculating coverage\n")
   MaxRanges <- GRanges()
   if(length(reads) > 0){
      AllCov <- coverage(reads)
      AllCov <- runmean(AllCov[names(AllCov) %in% ChrOfInterest],20)
      cat("Calculating Summits on ",ChrOfInterest," ..")
      Maxes <- which.max(Views(AllCov[[which(names(AllCov) %in% ChrOfInterest)]],ranges(bedRanges[seqnames(bedRanges) == ChrOfInterest])))
      if(any(is.na(Maxes))){ 
         NoSummitRanges <- bedRanges[is.na(Maxes)]
         Maxes[is.na(Maxes)]  <- (start((ranges(NoSummitRanges[seqnames(NoSummitRanges) == ChrOfInterest])))+end((ranges(NoSummitRanges[seqnames(NoSummitRanges) == ChrOfInterest]))))/2
      }
      MaxRanges <- GRanges(seqnames(bedRanges[seqnames(bedRanges) == ChrOfInterest]),IRanges(start=Maxes,end=Maxes),elementMetadata=elementMetadata(bedRanges[seqnames(bedRanges) == ChrOfInterest]))
      cat(".done\n")
   }
   return(MaxRanges)
}

getAnnotation = function(GeneAnnotation="hg19",AllChr){
   if(!is.null(GeneAnnotation)){
      if(GeneAnnotation == "hg19"){
         txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
      } else {
         stop('Unsupported annotation:',GeneAnnotation)
      }
      All5utrs <- reduce(unique(unlist(fiveUTRsByTranscript(txdb))))
      All3utrs <- reduce(unique(unlist(threeUTRsByTranscript(txdb))))
      Allcds <- reduce(unique(unlist(cdsBy(txdb,"tx"))))
      Allintrons <- reduce(unique(unlist(intronsByTranscript(txdb))))
      Alltranscripts <- reduce(unique(unlist(transcripts(txdb))))
      Promoters500 <-  reduce(flank(Alltranscripts,500,"start"))    
      Promoters2000to500 <-  reduce(flank(Promoters500,2000,"start"))
      LongPromoter20000to2000  <- reduce(flank(Promoters2000to500,20000,"start"))
      if(!missing(AllChr) & !is.null(AllChr)){
        All5utrs <- GetGRanges(All5utrs,AllChr=AllChr)
        All3utrs <- GetGRanges(All3utrs,AllChr=AllChr)
        Allcds <- GetGRanges(Allcds,AllChr=AllChr)
        Allintrons <- GetGRanges(Allintrons,AllChr=AllChr)
        Alltranscripts <- GetGRanges(Alltranscripts,AllChr=AllChr)
        Promoters500 <- GetGRanges(Promoters500,AllChr=AllChr)
        Promoters2000to500 <-  GetGRanges(Promoters2000to500,AllChr=AllChr)
        LongPromoter20000to2000  <- GetGRanges(LongPromoter20000to2000,AllChr=AllChr)        
      }
   }
   return(list(version=GeneAnnotation,All5utrs=All5utrs,All3utrs=All3utrs,Allcds=Allcds,Allintrons=Allintrons,Alltranscripts=Alltranscripts,
               Promoters500=Promoters500,Promoters2000to500=Promoters2000to500,LongPromoter20000to2000=LongPromoter20000to2000))
}



